Low-energy quantum dynamics of atoms at defects. Interstitial oxygen in silicon 
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The problem of the low-energy highly-anharmonic quantum dynamics of isolated impurities in solids 
is addressed by using path-integral Monte Carlo simulations. Interstitial oxygen in silicon is studied 
as a prototypical example showing such a behavior. The assignment of a "geometry" to the defect 
is discussed. Depending on the potential (or on the impurity mass), there is a "classical" regime, 
where the maximum probability-density for the oxygen nucleus is at the potential minimum. There 
is another regime, associated to highly anharmonic potentials, where this is not the case. Both 
regimes are separated by a sharp transition. Also, the decoupling of the many-nuclei problem into 
a one-body Hamiltonian to describe the low-energy dynamics is studied. The adiabatic potential 
obtained from the relaxation of all the other degrees of freedom at each value of the coordinate 
associated to the low-energy motion, gives the best approximation to the full many-nuclei problem. 
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I. INTRODUCTION 

The dynamics of impurities or, more generally, atoms 
at point defects in crystalline solids quite often gives rise 
to localized low-energy excitations, typically in the far 
infrared (FIR) spectral region, displaying patterns that 
reflect, in one way or another, a substantial deviation 
from the situation of atomic nuclei harnionically vibrat- 
ing around their potential minimum.EJ^H For example, 
interstitial H or Li around other substitutional impuri- 
ties in Si and Ge are delocalized among several position|S 
equivalent by symmetry, giving rise to FIR excitations 
Similar FIR excitation patterns are found in glasses 
Mixed crystals have recently been shown to exhibit a 
similar behavior.l3 

The substantial deviation from harmonicity in these 
systems, considering the quantum character of their dy- 
namics, is of great importance in their characterization. 
On one hand, the quantum character prevents the as- 
signment of a definite geometry to the defect structure, 
the nuclei showing high probability for being found away 
from the potential minimum location. On the other 
hand, the anharmonicity makes that the many-nuclei 
problem cannot be easily decoupled as usually in terms of 
normal modes of vibration, leaving open the question of 
what is the best decoupling approximation, that accounts 
for the experimental observations of those low-energy ex- 
citations. 

Interstitial oxygen (jQ^l.in silicon also displays nontriv- 
ial quantum behavior .BQu The oxygen atom is known to 
break a Si-Si bond, establishing two Si-0 bonds in the 
form of an oxygen bridge. The motion of the oxygen 
atom and its neighbors in the direction of the axis of the 
original Si-Si bond can be described in terms of harmonic 
vibrations, which account for some of the main features 
in the infrared absorption spectrum of the center .B It is 
the motion of the oxygen atom in the plane perpendicular 



to that axis what gives rise to the non-triyial behavior, 
characterized by a peculiar FIR spectrum.EJ The symme- 
try of the center corresponds to the D^d point group. 
This symmetry facilitates the a priori partial decoupling 
of the relevant degrees of freedom. It is this fact, together 
with the availability of experimental and theoretical in- 
formation, what makes this system particularly suitable 
for analysis of the questions raised above: (i) can we de- 
fine a geometry for the center, and (ii) what is the best 
way of decoupling the low-energy dynamics from the rest. 

The relevance of the geometry- question has already 
been pointed out in the literature.|j It was observed jLhat, 
for the effective potential obtained from experiment J3 the 
oxygen atom has maximum probability-density at the 
bond-center (BC) site in spite of being the potential a lo- 
cal maximum at that point. The arguments used there, 
however, are for the effective one-particle Hamiltonian, 
and have to be checked for the fully interacting problem. 

The importance of the decoupling problem is clear by 
the analysis of the experimental FIR spectra in terms 
of finding an effective one-particle potential in two di- 
mensions, such that an oxygen atom moving under that 
potential would reproduce the vibrational frequencies ob- 
served in the spectra. It is an axially-symmetric potential 
well with a local maximum at the bond-center and a min- 
imum around it, at 0.22 A off the BC site (mexican-hat 
shape) u The energy barrier amounts to ~ 1 meV. If we 
consider the fact that the oxygen atom is strongly in- 
teracting with its silicon neighbors, the question arises: 
what is the meaning of that potential? The multidimen- 
sional potential associated to all the nuclear degrees of 
freedom of the system is well defined, but not a two- 
dimensional one, unless a decoupling prescription is avail- 
able. In similar situations, when the harmonic decoupling 
is meaningless, propositions for decoupling are found in 
the literature ranging between two extremes. Arguing on 
the basis of the lightness of the impurity atom, some 
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authorsB propose that the effective one-particle potential for the motion of that impurity is the one obtained when the 
host (heavier) atoms are kept fixed at some defined positions (fixed-lattice potential; FLP, hereafter): the effective 
slower motion of the host atoms is presumed from mass considerations. The other extreme can be also stated in 
terms of the different times scales of the relative motions, but taking now into consideration the energies associated 
to those motions. In this latter case it is the impurity (or some of its degrees of freedom) the slow moving, in spite 
of its lightness. This argument leads to an adiabatic potential (AP) that is obtained by allowing the host atoms to 
relax for each value of the coordinates associated to the low-energy degrees of freedom. 

These questions are addressed in this paper by solving the full quantum many-nuclei problem of interstitial oxygen 
in silicon by means of path-integral (PI) Monte Carlo (MC) simulations. This paper is divided up in the following 
manner. In Sec. II we briefly describe the computational method and summarize the details concerning the potential 
used. Sec. Ill is devoted to the results and discussion. The geometry of the center Oi is assessed with the probability- 
density function of the oxygen nucleus. The decoupling of the low-energy motion is studied by comparing the results 
of the complete MC simulations with those derived according to different decoupling criteria. The paper closes with 
a summary (Sec. IV). 



II. COMPUTATIONAL METHOD 



A. Path-integral simulations 



The PI MC method has become a standard tool to study finite-temperature properties of quantum systems. In this 
Section we briefly present the details pJesant for the presentation of our results. More complete descriptions of the 
PI formalism can be found elsewhere. tHnlj The canonical partition function of P silicon nuclei plus the Q-j|mpurity 
can be expressed with the Trotter formula and the high-temperature approximation for the density matrixll3 as: 




The index j represents the path coordinate along a cyclic path, which has been decomposed into a number N of 
discrete intervals (Trotter number). Rj is a vector in a 3{P + I)-dimensional space, whose components are the 
Cartesian coordinates of the P + 1 nuclei (ri ...;rp+i,j). The cyclic condition for the path coordinate of each 
nucleus is expressed as R^v-i-i — Ri- The masses of the host and impurity atoms are mgi and toq , respectively; 
/3 = (fesT)"^, and fc^ is the Boltzmann constant. Eq. (^ coincides with the canonical partition function of a classical 
system with the effective interaction potential: 

TV 

i;off(Ri,...,RTv) =^{A(Rj,Rj+i)+iV-V(R,)}, (2) 

3 = i 



where 

p 



A(Rj, Rj+i) = I mo(rp+ij+i - rp+i^^f + ^ msi(rpj+i - Tpjf^ . (3) 



The index p indicates the particle, and goes from 1 to P for the silicon atoms, and takes the value P + I for the 
impurity. The function Weff(Ri, ■■■,'Rn) is the interaction potential of a classical system composed of P + 1 cyclic 
chains [one per nucleus; a total of N{P + 1) classical particles] characterized by a harmonic intrachain coupling with 
a force constant k = mN / 0^h^ [first term on the right-hand side of Eq. (|); m = msi or mo, depending on the 
nucleus] . Interchain coupling [second term on the right-hand side of Eq. (||) is restricted to those particles with the 
same index j, and this interaction is equal to that corresponding to the quantum particles, ViTlj), but renormalized 
by a factor (inverse of the number of discrete points along the path coordinate). The approximated expression 
for Z [Eq. ([^)] becomes exact in the limit N ^ oo, and is valid for distinguishable particles. This assumption of 
distinguishable particles is reasonable for the statistics of Si nuclei, since exchange effects are negligible. 

Equilibrium properties of the quantum system can be derived by Mejttppolis Monte Carlo sampling of the multidi- 
mensional integral associated to the partition function given byEq. (0)JlZm A simulation run proceeds via successive 
MC steps (MCS). In each MCS, the nuclei coordinates rpj are updated according to two different kinds of sampling 
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schemes. The first one consists on trial moves of the 
individual coordinates Vpj . The trials are performed se- 
quentially for every path coordinate j and every nucleus 
p. The second type of sampling ] corresponds to trial 
moves of the center of gravity of the cyclic paths, that 
are carried out sequentially for every nucleus p in the 
simulation cell. The number of MCS's employed for sys- 
tem equilibration was of 5 x 10^, while the calculation of 
ensemble average properties was performed over 2 x 10^ 
MCS's. For the Si atoms we have employed the average 
isotope mass of this element (m — 28.086 amu). The 
number N of discretized points for the path coordinate 
was made temperature dependent, and was taken as the 
integral number closest to 2000/r, a condition that guar- 
antees convergence in rthe total energy within a relative 
error smaller than 1%J23 



B. Potential 

The Monte Carlo simulations have been performed on 
a 2 X 2 X 2 supercell of the Si face-centered cubic (fee) 
cell containing 64 Si atoms and an oxygen impurity. The 
simulation cell was subject to periodic boundary condi- 
tions. The interaction between silicon atoms has been 
described by the rtliree-body potential developed by Still- 
inger and Weber .113 The potential between oxygen and sil- 
icon atoms has been designed to reproduce qualitatively 
the main features of the Oi defect: (i) the overall geom- 
etry, O breaking a bond between twa-,Si atoms; (ii) the 
observed low-energy (FIR) excitationsEI (even though the 
PI does not provide excited states, it does give internal 
energy versus temperature, which, at low temperatures is 
mainly controlled by these low-energy excitations); (iii) 
the vibrations of the cepljer at higher frequencies, known 
by infrared absorption;liil and (iv) the main features of 
the potential obtained from first-principles calculations- 
including bond lengths, bond angle, and Si relaxations.H 
This potential for the interaction between O and Si atoms 
is a function of the coordinates of both the impurity (ro), 
and the Si atoms (ri and r2) coordinated to the oxygen: 

y(ri,r2,ro) = K(rio) + T4(r2o) + V;(a) + V/(a) , (4) 

where r^o = fj — ro (* =lj 2), and a is the angle Si-O-Si. 
The potential functions are: 

Vr{r) = ^k{r - ref , (5) 
Vg (a) = si sin^ a + S2 sin^ a , (6) 

V;(q!) = ^i(c0Sa — COSOfe)^ , (7) 

with the following values of the constants: k = 35.6 eV 
A-2, re = 1.629 A, si = 1.49 eV, S2 = -0.7484 eV, 
li — 5.4 eV, and ag = 168°. The local geometry of the 



defect complex is shown in Fig.yJ. The absolute potential 
minimum corresponds to a geometry with the oxygen nu- 
cleus located at the off-center position M, at a distance 
of about 0.29 A from the Si-Si axis. The corresponding 
Si-0 distance is 1.52 A and the Si-O-Si angle is 158°. 
The relaxation of the nearest Si atoms along the [111] 
crystal direction amounts to 0.32 A. For comparison, the 
Si-0 distance obtained from total-energy Hartree-Fock 
calculationda of the Oi defect is 1.56 A, with an outwards 
relaxation of the nearest Si atoms of 0.36 A each. The vi- 
brational frequency of the Aig mode at the Oi center de- 
rived from our model potential is 587 cm~^, while the frfe 
quency reported in a cluster-Bethe-lattice investigational 
of the Oi defect was 569 cm^^. This mode corresponds 
to atom displacements along the Si-Si axis and has no 
infrared activity because of symmetry. The potential en- 
ergy for the impurity located at the BC site is higher 
than that found for the absolute minimum. This value 
agrees with that of ~ 1 meV corresgonding to the model 
potential of Yamada-Kaneta et al.iJ which was designed 
to reproduce the spectrum of low-energy excitations of 
the Oi center. When the oxygen is located at the BC 
site, our parametrized potential gives a relaxation of the 
nearest Si atoms of 0.35 A. 

In Fig. ^(a), a calculated potential energy surface is 
presented as a function of both the distance between the 
Si atoms nearest neighbors of O, rf(Si-Si), and the dis- 
tance from the oxygen atom to the bond-center site as 
oxygen moves in the plane perpendicular to the Si-Si 
axis, d(O-BC). For every point in the figure, the posi- 
tions of the other Si atoms were relaxed. The curves ob- 
tained by sectioning the surface for fixed values of d{Si- 
Si) correspond to FLP's. For distances rf(Si-Si) larger 
than 3.04 A the minimum energy is found when the oxy- 
gen atom is located at the BC site. However, at smaller 
values of (i(Si-Si) , the minimum is found for an off-center 
position. In Fig. ^(b) we present two sections of the en- 
ergy surface calculated at representative distances d(Si- 
Si). The broken line shows the FLP derived with the 
Si atoms fixed at their relaxed positions for the absolute 
minimum of the potential, the oxygen impurity being lo- 
cated at the off-center position M, with d(Si-Si) = 3 A 
(called FLP-1, hereafter). The energy barrier amounts to 
~ 24 meV. The solid line corresponds to a fixed-lattice 
potential (FLP-2) with the hosts atoms fixed at their re- 
laxed positions for the O nucleus located at the BC site 
and (i(Sii-Si2) = 3.05 A. The dotted line connecting the 
minima of both curves represents the adiabatic potential 
(AP) characterized by an energy barrier. 



III. RESULTS AND DISCUSSION 

A. Geometry 

In order to define a geometry for the Oi defect, we have 
studied the probability-density p{r) of finding the oxygen 
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nucleus at a distance r from the Si-Si axis along any di- 
rection by PI MC simulations at 10 K. The full quantum 
treatment included the O nucleus and the nearest and 
next-nearest Si neighbors (i.e., a total of nine atoms), 
while the remaining Si atoms in the simulation cell were 
fixed at their relaxed positions obtained for the system at 
the absolute potential minimum. The probability-density 
p{r) was normalized so that: 

oo 

J dr 2TTrp{r) = 1. (8) 



In Fig. ^ the density p{r) for the oxygen nucleus is shown 
by a full line. The maximum probability-density is found 
at the BC site, i.e., it does not correspond to the absolute 
minimum of the potential energy at the off-center site M. 
This peculiar behavior was already observed in a previous 
investigation assjjming an effective one-particle potential 
for the impurity.B The relevance of our calculation is that 
this nontrivial quantum delocalization of oxygen is con- 
firmed by a quantum approach for the full many-body po- 
tential. One expects that an impurity with larger mass, 
interacting through the same potential as that employed 
for oxygen, will approach a "classical" behavior, in the 
sense that those spatial regions with larger probability- 
density will correspond to regions of lower potential en- 
ergy. In Fig. ^, the function p{r) obtained by setting in 
the PI simulation the impurity mass toq — 60 amu, is 
displayed by a broken line. The maximum probability- 
density is found in this case at an off-center position. We 
have performed a series of simulations at 10 K varying 
the impurity mass in the range mo = 16-150 amu. In 
Fig. ^(b), the position of the maximum of the p{r) curves 
is presented as a function of the impurity mass mo . For 
a mass larger than about 50 amu we find a "classical" 
regime, where the position of the maximum of p{r) is 
off-center. 

We note that the assignment of a definite geometry 
to the Oi defect is difficult and may be meaningless, be- 
cause the probability-density function is very broad in 
the plane perpendicular to the Si-Si axis. This large 
spatial uncertainty is due to the zero-point motion of 
the impurity and therefore it cannot be reduced by de- 
creasing the temperature. From the point of view of a 
structure with minimum potential energy, which is the 
one usually adopted in total-energy investigations using 
electronic structure methods, the defect geometry would 
correspond to a non-linear arrangement of the oxygen 
and the nearest-neighbor Si atoms. This is of no physi- 
cal significance, however, since the maximum probability- 
density (and the symmetry) is for a linear disposition of 
the atoms. It is the highly anharmonic situation found 
for oxygen in silicon that leads to contradictory pictures 
concerning the defect geometry. This point could be rele- 
vant also for some crystalline phases, like /3-cristobalite, 
where different structural models have been proposed, 
and some controversy has arisen concerniM-the existence 
of linear Si-O-Si units in this struct ure.EilE3 



B. Decoupling into a one-particle problem 

The results derived from the full quantum mechanical 
treatment of the Oi defect can be used to test the quality 
of several alternatives (FLP versus AP) for the decou- 
pling of the many-nuclei degrees-of-freedom into a one- 
body potential. In Fig. |^(a), the "exact" probability- 
density p(r) for the oxygen nucleus obtained by the full 
PI simulation at 10 K is represented by open squares 
and is compared to the curves derived by four different 
decoupling schemes. We have analyzed three different 
fixed-lattice potentials FLP and the adiabatic potential 
AP. For each one of these potentials we have performed 
a one-particle PI simulation to obtain p{r) at 10 K. The 
broken line in Fig. ^a) is the result derived from the 
potential FLP-2, defined with the Si atoms fixed at the 
relaxed positions obtained with the O impurity located at 
the BC site. The maximum probability-density is found 
at the BC site in agreement with the result obtained for 
the full potential. However, the probability-density away 
from the BC site decays too fast and the curvature of the 
p{r) curve does not compare well with the result of the 
full PI simulation. The dashed-dotted line corresponds 
to the potential FLP-1, defined with the Si atoms fixed 
at the relaxed positions obtained with the O impurity 
located at the off-center site M. This potential FLP-1 
leads to an even more unrealistic probability-density, as 
the maximum density is found close to the off-center po- 
sition M. The dotted line in Fig. ^(a) was derived from 
the third choice of a FLP, where the Si atoms were fixed 
at their equilibrium positions obtained from the MC tra- 
jectory generated by the full quantum simulation of the O 
impurity and the Si atoms. This potential FLP-3 gives 
a better agreement with the full simulation. However, 
in spite of its much larger computational requirements 
(it needs the results of the full quantum treatment), the 
improvement is still unsatisfactory. This and the other 
FLP approaches discussed above are not able to repro- 
duce closely the density distribution of the impurity. At 
last, we present in Fig. ^(a) the results found for the oxy- 
gen impurity moving in the adiabatic potential defined in 
the plane perpendicular to the Si-Si axis [dotted line in 
Fig. 1(b)]. This pir) curve is shown by a full line and 
agrees closely with that derived by the full quantum sim- 
ulation. We conclude that the best decoupling scheme 
for the Oi defect, in the plane perpendicular to the Si-Si 
axis, corresponds to the AP approximation. 

It is interesting to test whether the AP potential is also 
able to reproduce the results of the full simulation shown 
in Fig. H(b), by a series of one-body simulations for dif- 
ferent impurity masses. The results obtained in these PI 
MC simulations for the position of the maximum of p{r) 
are given in Fig. ^(b). The agreement with the data of 
the full quantum simulations of Fig. ^(b) is good. In par- 
ticular, the crossover between impurity- mass ranges with 
maximum probability-density at BC or off-center, lies at 
about 50 amu, as found in the full quantum problem. 
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As the mass of the impurity increases, its motion with 
respect to the host atoms will become relatively slower, 
therefore the adiabatic potential remains a good approx- 
imation to the many-body problem. 

Finally, we compare in Tab. | some spectroscopic 
information, obtained by solving numerically the two- 
dimensional Schrodinger equation of an oxygen nucleus 
moving in the AP_and FLP-1 potentials, with available 
experimental data.Ej The frequencies derived from the AP 
potential show a closer agreement to the experimental re- 
sults than the frequencies calculated with the FLP-1 po- 
tential. The AP transitions are about 40 % larger than 
those derived from experiment. However, the renormal- 
ization of the transitions to the value of the first exci- 
tation reduces this discrepancy to 12 %. This deviation 
from experiment can be due either to limitations in the 
adiabatic approximation or to shortcomings in the em- 
ployed potential model. The main conclusion from the 
data in Tab. | is that the adiabatic potential appears 
again as a better decoupling approach to the full many- 
nuclei problem than the fixed lattice potential. 



IV. CONCLUSIONS 

The definition of a defect geometry for an isolated oxy- 
gen impurity in silicon is conditioned by a substantial de- 
viation from harmonicity, as the configuration with min- 
imum potential energy does not correspond to a maxi- 
mum in the probability-density function of the oxygen 
nucleus. The minimum energy configuration is found for 
a bent Si-O-Si configuration. However, the probability- 
density function of the impurity, as derived by PI MC 
simulations, displays the symmetry of a linear arrange- 
ment of atoms. This is a case in which the structure 
associated to the minimum energy is not the best choice 
to derive some physical properties of the defect complex. 
The anharmonicity of the zero-point motion of the impu- 
rity and the lattice atoms reduces the importance of the 
minimum-energy structure for defining the properties of 
the impurity center. 

An important property of the Oi center is that the 
dynamics in the plane perpendicular to the Si-Si axis is 
characterized by frequencies of about 30 cm^^. These 
extremely low frequencies, even for a relatively light im- 
purity as oxygen, cause that the best decoupling scheme 
for the reduction of the full many-body problem into a 
one-body treatment is given by the adiabatic potential, 
where for each impurity position the Si atoms (moving 
faster) relax to their equilibrium sites. 

Low-energy excitations are normally associated to lo- 
calized defect states. Some of the physics, however, is 
very much related to other dynamical phenomena like 
diffusion. Specifically, the problem of the delocalization 
of light impurity atonas among symmetry-equivalent sites 
around another atomQ is essentially the same as what is 
found for the diffusion of those light atoms, where the 



delocalization is among sites equivalent by translational 
symmetry. The discussion in this paper about the defi- 
nition of a static defect structure is partially transferable 
to the dynamical problem of defining diffusion paths for 
quantum particles. 
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FIG. 1. Schematic representation of the atom disposition 
in the Oi center. Alternative sites for the oxygen nucleus are 
indicated by BC (bond-center) and M (off-center). The abso- 
lute minimum of the employed potential model corresponds 
to the off-center site M, 0.29 A away from the BC site. 



TABLE L Observed and calculated FIR transitions 
(cm^^) of oxygen in silicon. The calculated values correspond 
to the AP and FLP-1 potentials. In parenthesis are given the 
relative values with respect to the first transition. The ob- 
served values are taken from Ref. [1]. 



Transition 


obs. 


AP 


FLP-1 


|0,0)^|1,±1) 


29.3 (1) 


37.1 (1) 


24.1 (1) 


|1,±1)^|2,±2) 


37.8 (1.3) 


51.6 (1.4) 


51.6 (2.1) 


|2,±2) |3, ±3) 


43.3 (1.5) 


62.9 (1.7) 


70.9 (2.9) 


|1,±1)^|2,0) 


49.0 (1.7) 


72.6 (2.0) 


140.3 (5.8) 



FIG. 2. (a) Potential energy for the nuclear motion shown 
versus the O distance to the BC site within the plane per- 
pendicular to the Si-Si axis, and versus the Si-Si distance 
between host atoms nearest- neighbors of O. The contour in- 
terval is 3 meV. (b) Three different one-dimensional cuts are 
presented: (i) for the Si atoms fixed at the positions of the ab- 
solute potential- minimum (fixed- lattice potential FLP-1, bro- 
ken line); (ii) for the Si atoms fixed at their relaxed positions 
for O at the BC site (FLP-2, full line) , and (iii) for the Si 
atoms relaxing at each O position (adiabatic potential, AP, 
dotted line), (i) and (ii) are cuts of the potential surface in 
(a), in vertical planes parallel to the d(O-BC) axis. 



FIG. 3. (a) Probability-density p(r) for the O nucleus ob- 
tained in the full PI simulation at 10 K. r is the distance from 
the impurity to the Si-Si axis. The full line corresponds to 
mo = 16 amu and the broken line to mo = 60 amu. (b) 
Distance from the probability-density maximum to the Si-Si 
axis, as a function of the impurity mass. The broken line is a 
guide to the eye. 



FIG. 4. (a) Comparison of probability-density curves of the 

O nucleus for different one-particle potentials with that de- 
rived by the many-particle simulation (open squares). The 
broken, dotted, and dashed-dotted lines correspond to differ- 
ent FLP models (see text), while the full line is the result of 
the adiabatic potential (AP). (b) Position of the maximum 
in the probability-density p(r) corresponding to the adiabatic 
potential, as a function of the impurity mass. The broken line 
is a guide to the eye. 
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